
NASA Contractor Report 4270 








NASA Contractor Report 4270 


A Direct-Inverse Method 
for Transonic and Separated 
Flows About Airfoils 


Leland A. Carlson 
Aerospace Engineering Department 
Texas A&M University 
College Station , Texas 


Prepared for 

Langley Research Center 

under Grant NSG-1174 


I\J/\SA 

National Aeronautics and 
Space Administration 

Office of Management 

Scientific and Technical 
Information Division 


1990 



A DIRECT -INVERSE METHOD FOR TRANSONIC AND SEPARATED 
FLOWS ABOUT AIRFOILS 


Leland A. Carlson 
Aerospace Engineering Department 
Texas A&M University 


SUMMARY 


A direct- inverse technique and computer program called TAMSEP that can 
be used for the analysis of the flow about airfoils at subsonic and low 
transonic freestream velocities is presented. The method is based upon a 
direct- inverse nonconservative full potential inviscid method, a Thwaites 
laminar boundary layer technique, and the Barnwell turbulent momentum integral 
scheme ; and it is formulated using Cartesian coordinates. Since the method 
utilizes inverse boundary conditions in regions of separated flow, it is 
suitable for predicting the flowfield about airfoils having trailing edge 
separated flow under high lift conditions. Comparisons with experimental data 
indicate that the method should be a useful tool for applied aerodynamic 
analyses . 
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isentropic speed of sound 

boundary layer coefficient in separated pressure correlation 

airfoil drag coefficient 

airfoil lift coefficient 

pressure coefficient 

Mach number 

velocity 

velocity components in the x and y directions, respectively 
transformed velocity at boundary layer edge 
velocity in the boundary layer 

law- of - the -wall and law- of - the -wake velocity parameters 
Cartesian coordinates 
angle of attack 

ratio of specific heats, assumed to be 1.4 
circulation 

boundary layer thickness 
polar coordinate 
computational coordinates 
potential function 
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perturbation potential 


Subscripts : 


00 

freestream condition 

b 

body 

e 

boundary layer edge 

1 -j 

grid location 

LE 

leading edge 

SEP 

separation point 

TE 

trailing edge 

if n 

differentiation 

y 

differentiation 


GENERAL DESCRIPTION 


INTRODUCTION 

1-3 

Over the past decade, several finite -difference potential flow methods 

have been developed and successfully used for the design and analysis of 

subsonic and transonic airfoils at and near cruise conditions. However, in 

the analysis of high performance airfoils, aerodynamic is ts would also like to 

be able to predict airfoil pressure distributions and aerodynamic coefficients 

at high lift, high angle-of -attack conditions. Since such situations are 

frequently characterized by regions of separated flow on the upper surface and 

are dominated by strong viscous interaction effects, inviscid methods alone 

are not applicable. Furthermore, subsonic- transonic analysis methods^'^ which 

couple inviscid and boundary layer solutions typically only include the 

effects of weak viscous interaction and generally fail to give accurate 

results when separated flow exists on the upper surface. 

5-8 

However, it has been demonstrated that the direct- inverse technique 
coupled to a suitable boundary layer method can be successfully applied to low 
speed flows about airfoils having massive separation. In addition, Barnwell^, 
Dvorak and Choi^, and Taverna^ have developed similar methods for transonic 
flows. Barnwell's method, however, Is limited in application in that it 
utilizes for its inviscid solver the transonic small perturbation equation. 
Further, references 9 and 11 only include the effects of viscous Interaction 
due to a turbulent boundary layer. On most airfoils, particularly on the 
lower surface at high angles of attack, extensive regions of laminar flow 
exist . 

This report describes a flow model and computer program, called TAMSEF, 
which can be used to predict the flowfield about a single element transonic 
airfoil at high angle of attack high lift conditions with trailing edge 
separation. Since the method is based upon the TRANDES^ and TRANSEP^ codes, it 
can also be used for subsonic- transonic analyses not involving separation. 

METHOD OF APPROACH 

The present approach is based upon the direct- inverse method developed in 
the TRANDES and TRANSEP programs and the ability of this method to use either 
the displacement surface (airfoil ordinate plus displacement thickness) or 
pressure as the airfoil boundary condition. For the high angle-of -attack 
case, the airfoil lower surface only experiences weak viscous interaction and 
frequently has a long laminar run before transitioning to fully turbulent 
flow. Thus, the present model includes an initial laminar boundary layer 
calculation in its viscous interaction section. On the upper surface the 
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boundary layer is also initially laminar, but it quickly becomes turbulent in 
character followed in many cases by boundary layer separation and a separated 
zone which can extend over a significant portion of the airfoil surface. In 
the present model, this separated region is treated inversely in that the 
pressure distribution along the effective displacement surface streamline is 
determined iteratively as part of the solution and used as the airfoil 
boundary condition. Consequently, the present method has been modeled as 
shown on figure 1. 

To obtain the inviscid portion of the flowfield, the full potential 
equation for two-dimensional compressible flow is used in nonconservative form 

as 

(« 2 - *x)*xx ■ 2 W*y + ( a2 • *y)*yy " 0 (1> 

where the subscripts denote partial differentiation. By defining a 
perturbation potential, 4> < such that 

$ = xq^cos a + yq^ sin a +q a <f> (2) 

where the velocity components are given by 

U - 3> x - q^Ccos a + </> x ) 


and 


v “ $ y - q^sin a + <f> y ) 

the governing equation in terms of the perturbation potential 
as 


(3b) 

can be written 



- u2 Kx • 2UV *xy + ( a2 - 


v2 )*yy - 0 


(4) 


with 


“ 2 - a 2 - N" 1 ) ( U2 + v2 ' q -) 

The nonconservative form of the potential equation was selected for the 
present problem because results for two-dimensional flows obtained with it 
agree better with Euler solutions than those obtained using the fully 
conservative form of the equation^ . In addition, the conservative 
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formulation appears to break down in two-dimensional cases shortly after the 

13 

onset of supercritical flow 

In the present model, equations (3-5) are finite differenced using a 
rotated difference scheme and solved iteratively using column relaxation in a 
stretched Cartesian grid which maps the infinite domain to a finite 
computational box. The appropriate boundary condition at infinity Is^'^ 


1 - (8 - «> ( 6 ) 


S tan ' 1 [ 


where 6 is the polar angle, and T is the circulation, which Is determined by 
the change in potential across the Kutta- Joukowski cut at the trailing edge of 
the airfoil. 

Likewise, the appropriate airfoil boundary condition in the direct 
regions (regions without separation having only weak viscous interaction) is 
the flow tangency condition given by the ordinates of the airfoil displacement 
surface, i.e. 


dy 

dx 



q«l 

r sin a + ^yb 

<loO 

cos a + <£ xb 


(7) 


In the inverse or separated flow region the pressure distribution along the 
effective displacement surface streamline is considered specified and used as 
the boundary condition. As shown in reference 2, this approach leads to a 
derivative boundary condition for the inverse region of the form 


<f, xh = -cos a + 


1 + 


1 - 


1 + 


7m:c 




,i 

2 


°Q Pbl 


(t - dm: 


( 8 ) 


Complete details concerning the finite difference scheme, the stretched 
Cartesian grid system, and the treatment of the boundary conditions are given 
in references 1, 2, and 4. 

To include viscous effects, the basic approach is to calculate a boundary 
layer displacement thickness for the weak interaction regions and to use It to 
correct the location of the displacement surface (i.e., airfoil ordinate plus 
displacement thickness) . For the strongly interacting separated zone on the 
upper surface, the pressure is determined from the interaction solution and 
the location of the displacement surface Is computed by integrating the 
surface tangency condition, equation (7) , with the initial conditions 
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specified by the displacement surface ordinates at the separation point, which 

is the interface between the two regions. The location and slopes of the 

displacement surfaces are updated regularly throughout the iterative solution. 

In the present method, the laminar portion of the boundary layer is 

computed using a compressible Thwaites method similar to that used previously 

in TRANSEP 6 . The transition location is determined from a Granville type 

correlation 16 based upon the difference between the local momentum thickness 

Reynolds number and the value at the laminar instability point combined with 

the pressure gradient history. Sometimes, particularly on the upper surface 

at high angles of attack, laminar separation is predicted upstream of the 

transition point. In these cases, the local momentum thickness Reynolds 

number Is compared to an empirical correlation in order to determine if the 

laminar bubble is long or short. If the bubble is short, its length is 

assumed to be one horizontal delta-x grid width and the turbulent flow 

computation is initiated at the next downstream grid point. If the estimate 

indicates that the bubble is long, the calculation proceeds, but a warning is 

printed which indicates that the results are probably in error. 

After transition, the turbulent boundary layer is computed using the 

simplified Kuhn and Nielson method (SKAN) as developed by Barnwell In 

reference 9. This method was selected because it is efficient, reliable, and 

yields excellent predictions of displacement thicknesses and separation point 

location. The SKAN turbulent boundary layer method solves the integral forms 

of the momentum equation, moment of momentum equation, and the derivative of 

the Coles' law-of- the -wall law-of-the wake relationship applied at the 
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boundary layer edge. After considerable effort , these equations can be 
transformed into a set of simultaneous ordinary differential equations, I.e. 

"A 1 

du du* d<$ dU e du e a^ 

a il dx~ + a i2 dx + a i3 dx “ b i du^ dx~ + c i 

At 

which can be solved for the wall friction velocity, u , the wake parameter u^, 
and the boundary layer thickness 6, using a second-order predictor-corrector 
technique. The remaining quantities of interest such as displacement 
thickness and momentum thickness are then determined from these variables. 

The numerical integration Is terminated at the separation point, where the 
wall friction velocity, u , vanishes. 

The method uses a two -layer eddy-vi scosity model which utilizes the 
Prandtl mixing length concept in the inner region and a Clauser model in the 
outer zone. In addition, the method assumes an adiabatic wall, ignores the 
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laminar sublayer, and approximates the intermittency factor as well as several 
density ratios appearing in the fundamental equations. 

Thus, on the lower surface the flow is computed using direct boundary 
conditions (airfoil specified) including the effects of weak viscous 
interaction. On the upper surface, the flowfield is also computed directly 
with viscous interaction up to the separation point, which is determined as 
part of the boundary layer solution. Downstream of separation, inverse 
boundary conditions are utilized, and the pressure must be specified. 
Fortunately, if the skin friction at the wall is assumed to be zero in the 
separated zone, the SKAN formulation can be used to obtain a closed form 
solution for the velocity, and hence the pressure, at the outer edge of the 
separated zone. The resultant analytic expressions for the velocity and 
pressure are 



( 10 ) 

(ID 


As can be seen, the separated pressure depends upon the flowfield 
solution via the inviscid perturbation potentials at the separation point and 
the trailing edge, the size of the separated zone, and (through c^) the 
boundary layer solution at the separation point. In addition, this closed 
form solution predicts a variable pressure distribution for the separated 
region. At low freestream Mach numbers this variation is extremely small and 
is essentially constant. However, at freestream Mach numbers of 0.3 and 
above, the variation becomes significant and influences the resultant 
flowfield solution. This trend and the separated pressure variation is in 
accord with experimental observations and is a significant improvement over 
previous methods which assumed constant pressure in the separated zone 
regardless of the flow conditions. At low speeds the separated pressure is 
essentially constant and the complexity introduced by equations (10) -(11) may 
not be warranted. Thus, the present method contains the option of either 
using a constant pressure in the separated region or the variable distribution 
determined by the closed form solution given above. 

In principle, the separated region and the wake should be accurately 
modeled with respect to physical phenomena and Internal details. 

Consequently, several other investigators ^ ’ ^-18 ^ ave attempted to develop 
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detailed descriptions of these regions. In the present model, however, the 
wake region contains very few computational points because the coordinate 
system rapidly stretches to infinity. Thus, the wake is assumed to be 
inviscid with a constant pressure trailing edge formed by the upper and lower 
displacement surfaces. Fortunately, extensive numerical experiments with the 
present and previous^ models Indicate that the pressure distribution and 
aerodynamic coefficients are primarily dependent upon obtaining accurate 
predictions for the location of the separation point and the magnitude and 
variation of the separated pressure. Apparently, for many problems of 
interest the details of the wake region are of secondary importance. Since 
the present method obtains the separation point location directly from the 
solution for the wall friction velocity, u , and the pressure variation from a 
solution which couples the inviscid and viscous parts, it should yield 
reasonable engineering results. 

INTERACTION AND ITERATION PROCEDURE 

The iteration and interaction procedure used in TAMSEP is similar to that 
used In the low- speed program TRANSEP and is outlined in schematic form on 
figure 2. The program first reads all necessary Input and initializes the 
perturbation potential at all grid points to zero. Next it computes 
transformation factors and coordinates associated with the stretched Cartesian 
grid for the initial grid specified by the Input data. Included in this 
process is the computation of all airfoil ordinates and slopes required on the 
computational grid. 

Since the initial grid is normally very coarse with a default size of 
13x7, only fifty inviscid iterative cycles are computed on this grid. The 
calculation procedure used for the inviscid potential equation is the same 
iterative successive column over-relaxation scheme used in TRANDES ’ . This 

limited number of cycles serves to rapidly create an approximate starting 
solution for succeeding grids. After these are completed, the grid is then 
halved and the solution interpolated onto a new grid, which has a default size 
of 25x13. 

After obtaining all necessary coordinates , stretching factors , airfoil 
ordinates and slopes, etc., for this second grid, the method then performs 
fifty inviscid iterative cycles before considering any type of viscous 
interaction. Experience has shown that it is important to perform a limited 
number of inviscid cycles at the beginning of each new grid in order to 
eliminate any "problems" introduced by grid halving and interpolation. 

After these initial iterations, the program then checks to see whether or 
not the user desires viscous interaction to be included by examining the value 
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of the variable ITACT. If viscous interaction is desired, which is specified 
by the ITACT default value of one, the program then checks to see if an 
initial laminar boundary layer Is to be included (ILAM=1) or if the viscous 
calculations are to be for a turbulent boundary layer with user specified 
transition points (ILAM=0). 

Upon completing the boundary layer computations for the current flowfield 
solution, the program then calculates the ordinates and slopes of the upper 
and lower displacement surfaces. Since it only involves weak viscous 
interaction, the lower surface computations are from the leading edge or the 
lower surface stagnation point, whichever is further aft, to the trailing 
edge. However, on the upper surface they are only from the leading edge to 
the separation point or to the trailing edge, whichever is less. This process 
involves smoothing of the displacement thickness values, properly adding them 
to the airfoil ordinates, and spline fitting the resulting points. 

At this point the procedure depends upon whether or not separation has 
been detected on the upper surface. If separation does not exist prior to the 
last grid point on the airfoil upper surface, additional inviscid cycles are 
performed before returning to the viscous interaction loop. However, if 
separation is predicted, then the method must determine the pressure 
distribution and the location of the displacement surface in the separated 
zone . 

Exactly how the separated pressure distribution Is determined depends 
upon the user specified variable KSEP. If KSEP is zero, the pressure is 
assumed to be constant in the separated zone and is computed by 

r ' 2 f^TE ' ^SEpI 

c p x — r~£ (12) 

SEP X TE X SEP 

While this expression is a small perturbation approximation for C p g Ep , its 
usage has been found to be accurate and adequate for low speed incompressible 
flows. At frees tream Mach numbers of 0.3 and higher, however, the variable 
separated pressure option, specified by a KSEP value of one, should be used. 

In this case, the pressure distribution along the displacement streamline in 
the separated region is determined by equations (10) and (11) above. Note 
that both approaches determine the separated zone pressure, which depends upon 
the current solution, by conditions at both the separation point and at the 
trailing edge and not just on conditions in the vicinity of separation. This 
result is in agreement with the conclusion of Gross 19 that conditions at the 
downstream end of the separation zone partially influence the separation 
pressure level. 
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After determining the separated pressure distribution corresponding to 
the predicted separation point and the current potential flow solution, the 
corresponding upper surface displacement surface must also be computed for the 
separated zone. When separation exists, the previous method of adding the 
computed displacement thickness to the original airfoil ordinates is 
inappropriate since the values for displacement thickness predicted by the 
SKAN method are probably inaccurate in separated regions. Instead, the 
present approach is to solve, using the current potential flow solution, the 
differential equation 


dy V 
dx b " U b 


sin a + g b * wb 
cos q + f b 4>£ b 


(13) 


for the y-ordinates of the separated displacement surface as a function of x 
Based upon previous studies * , equation (13) is solved using the Runge-Kutta 
method of order four and the displacement surface ordinate at the separation 
point as the initial condition. In addition, in the process of solving this 
equation <f > ^ and 4 > ^ must be evaluated by finite differences. While several 
formulations are possible^, numerical studies indicate that accurate 
displacement surfaces are obtained using the following 




(14a) 



(14b) 


In equations (14), the point (i,j-l) is the first ghost point below the 
displacement surface. Its value is determined as part of the inverse pressure 
boundary condition. 

Since the present process is iterative and the potential solution uses 
the separated pressure distribution as an inverse boundary condition, the 
solution of equations (13) and (14) should yield upon convergence a separated 
zone displacement surface or free streamline that is compatible with the 
pressure distribution and potential flow solution. 

At this point in the iteration- interaction procedure a check is made to 
see if the solution has converged or if the maximum number of iterations for a 
given grid size has been exceeded. If neither situation is true, ten more 
inviscid cycles with the new displacement surfaces and separated pressure 
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distribution are performed prior to repeating the viscous interaction loop. 

If, however, either condition is satisfied and the finest grid specified by 
the user has not been used, the grid is refined and the entire process shown 
on figure 2 is repeated. If the last grid solution has been obtained, then a 
final output is printed and the solution is finished. 

It should be noted that the calculations on a given grid are stopped and 

assumed to be converged when the maximum perturbation potential change is less 

than some user specified value. However, when separation is present it is 

usual for the calculations on each grid to be terminated due to the number of 

iterative cycles exceeding a maximum user specified value (particularly on 

computers which only retain seven significant digits). In those cases, the 

existence or degree of convergence can be determined by examining the 

variation in the number of supersonic points, the location of separation, and 

the trailing edge ordinate of the upper displacement surface. All these 

values are printed out every ten Iterative cycles. If they stabilize prior to 

the end of the computation on a given grid, then the results can be assumed to 

be converged. Normally, it is sufficient to perform 800 cycles on the coarse 

grid (25x13), 400 on the medium grid (49x25), and 400 on the fine grid 

(97x49), although occasionally more may be needed. In determining 

convergence, it should be remembered that the present method is supposed to 

obtain a steady state solution. At angles of attack above maximum lift, the 

20 21 

actual flowfield about an airfoil is usually unsteady * . In those cases, 

the present method probably will not converge and may enter some type of 
oscillatory behavior which appears to represent an unsteady flow pattern. 
However, the present method Is not "time -accurate" and such results should 
only be viewed as indicative of the presence of significant unsteady 
phenomena . 
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USER INSTRUCTIONS 


CODE DESCRIPTION 

The TAMSEP code consists of a main program and eighteen subroutines. The 
subroutine and their relationships are shown in a subroutine tree on figure 3. 
The subroutine names and their functions are as follows: 

FOIL Reads in initial airfoil shape and determines ordinates and slopes at 
computational points. 

VISACT Computes turbulent boundary layer when viscous interaction included. 

THWAIT Computes the laminar boundary layer when viscous interaction 
included. 

FIT2 Curve fit routine used by Thwaites method. 

VALUE Initializes the flowfield to zero perturbation potential. 

SOLVE Sets up the matrix coefficients used in the SLOR relaxation scheme. 
PRESS Computes the pressure distribution on the airfoil. 

COORD Sets up the coordinates in the computational and physical grids and 
computes the stretching factors. 

FL0W1 Solves flowfield in front of the airfoil. 

FLOW 2 Solves flowfield in the direct region above and below the airfoil. 

FL0W3 Solves flowfield in the inverse region. 

WAKE Solves the flowfield behind the airfoil. 

SHAPE Computes the shape of the airfoil displacement surface in the 
separated zone. 

TRID Tridiagonal equation solver. 
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HALVE 


Doubles the grid size and interpolates old values to obtain starting 
values on the new grid, 

PLOT Creates a printer plot of the Cp and displacement surface. 

ARC Determines the arc length of the airfoil coordinates and splines the 

coordinates versus arc length. 

SPLINE Computes a cubic spline through a set of points. 

BDLY Computes the turbulent boundary layer at the end of an Inverse design 

calculation. 

The TAMSEP code is written FORTRAN IV programming language and is 
designed for use on IBM, AMDAHL, CDC, DEC, and similar computers. In 
nonoverlay mode it requires less than 320,000 bytes on an Amdahl 470/V8. 

Using a FORTRAN H extended compiler at the optimization level two, it needs 
about 17 seconds for compilation and obtains a solution on a 97x49 grid in 
about 160 seconds at a rate of around 15,000 points/second. Some slight 
modification to formats, etc. may be required to run the program on different 
computer systems or under a FORTRAN 77 compiler. 

INPUT DESCRIPTION 

The Input to the TAMSEP code is read in eight separate blocks. The first 
one contains a user supplied title, while the second and third blocks specify 
all the floating point and integer parameters needed to run the program. 

These parameters are input via namelists and if not specified are assigned 
default values by the program. Blocks four and seven are optional and are 
only included when the parameter IREAD is one. They read a non-zero 
perturbation potential starting solution and an initial airfoil description. 
Blocks five and eight are associated with input for the design option in the 
program and are only included when the parameter INV is one. Finally, block 
six contains the description of the airfoil under consideration. For an 
analysis computation, only blocks one, two, three and six would be included In 
the input stream. 

DETAILED INPUT DESCRIPTION 

Input Block 1: Title. 

This block consists of a single line of input and is read by the main 
program. 
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NTITLE Description of case. Up to 80 alphanumeric characters. Appears on 
the printed output at the beginning of the results of each grid. 

Input Block 2: Floating Point Parameters, 

This block of input is read by the main program via a namelist called 

FINP. 

M Freestream Mach number (real variable). Default 0.5 

W Relaxation factor for subsonic points. Should be in the range 

0 < W < 2.0. Default 1.7 

XI X location where the direct mode calculation procedure stops. In the 

analysis mode it should be set to 0.5 (i.e. the trailing edge). In 
the inverse (design) mode it is usually set to slightly less than the 
third point from the leading edge or larger. Default 0.5 

X2 End of the inverse region. For analysis cases set to a large number. 

In the inverse (design) case set to 0.5 (i.e. the trailing edge). 
Default 10000.0 

ALP Angle of attack in degrees. Default 0.0 

EPS Subsonic damping factor to match difference equations at sonic line 

if needed. Has no effect on the accuracy of the solution. Only 
affects stability and convergence rate. Normally it is not needed. 
Default 0.0 

EPSS Supersonic damping factor for iterative stability. Has no effect on 

the accuracy of the converged solution, only on the stability and 
convergence rate. Should typically be about M^ ax - 1, where M max is 
the maximum local Mach number. Default 0.4 

X4 The positive X locations where the coordinate stretching changes. It 

should be near the airfoil trailing edge. Default 0.49 

S4 The positive psi value in the computational plane where the 

stretching changes. Default 2.0 
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CONV 

A1 

A2 

A3 

RN 

XIBDLY 

CIR 

CDCORR 


Convergence criteria control value. Iterations stop when the maximum 
change in the perturbation potential between relaxation cycles is 
less than CONV. Default l.E-05 

Stretching constant for the Y direction. It can be used to control 
the Y and eta spacing near the horizontal axis. It is usually best 
to have the psi and eta spacing equal near the leading edge of the 
airfoil. Default 0.246 

First stretching constant for the X-direction. It is equivalent to 
( 2 /tt) (dx/dO at f The value of A2 determines the horizontal 

step size near the leading and trailing edges, i.e. 

ay tt A2 A2 2(1 + S4) 

x=x 4 2A£ "2 (IMAX-1) 

See Appendix A of reference 4. Default 0.15 

Second stretching constant for the x-direction. It determines the 
physical location of the vertical grid line adjacent to the grid side 
edge. Default 3.87 

Frees tream Reynolds number based on chord length. Used only when 
viscous interaction (ITACT-1) included. Default 20E+06 

The x- location at which upper surface transition is assumed to occur. 
The turbulent boundary layer calculation starts at the next grid 
point. The relationship to percent chord is: 

XIBDLY « (%chord - 50.0)/100.0 

Used only if viscous interaction included (ITACT=1) and laminar 
boundary layer ignored (ILAM=0) . Default -0.44 

Circulation about airfoil. If an initial solution is input (IREAD-1) , 
it must be the corresponding value of circulation (CIR=CL/2 . 0) . 

Default 0.0 

Correction to the wave drag coefficient. Because of the lack of a 
large number of points In the leading and trailing edge regions, the 
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RDEL 

RDELFN 

SP 

XSEP 

RCPB 

CPB 

XMON 

XPC 


XLBDLY 


wave drag coefficient has an error associated with grid size, 
spacing, and lift coefficient. The magnitude of CDCORR must be 
determined by the user by empirical methods. Note that the 
correction should be different for each airfoil and grid combination. 
Default 0.0 


Relaxation parameter for the boundary layer displacement thickness. 

It is used only when viscous interaction is included (ITACT=1) and 
IMAX is less than or equal to 55. Default 0.25 

Fine grid relaxation parameter for the boundary layer displacement 
thickness. It is used only when viscous interaction is included and 
IMAX is greater than 55. Default 0.125 

Maximum value allowed for the Nash-Macdonald separation parameter 
when x < XSEP. Used only in the design case (INV=1) when computing 
the boundary layer over the design surface. Default 0.004 

X location after which the Nash -Macdonald separation parameter can 
exceed SP. Used only in the design case (INV=1) when computing the 
boundary layer over the design surface. Default 0.44 

Not used. Ignore. 

Not used. Ignore. 

Not used. Ignore. 

Location after which the lower surface displacement thickness is 
required to continue decreasing once it has started to decrease. 
Upstream of XPC the displacement thickness is required to be 
monotonically increasing. For most aft-cambered airfoils it should 
be set to 0.1, and for conventional airfoils it should be set to 0 . 5 
Default 0.1 

The x- location at which lower surface transition is assumed to occur 
Same relationship to chord as XIBDLY. Used only if viscous 
interaction included (ITACT^l) and laminar boundary layer ignored 
(ILAM=0) . Default -0.44 
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RLAX 


Relaxation parameter for the separated pressure level in the constant 
separated pressure option (KSEP=0) . Sometimes needed to enhance 
convergence. Used only when ITACT=1, IMASS=1, and KSEP-0. Default 
1.0 

RADUS Leading edge radius of the airfoil nondimens ionalized by the chord. 
Used only if ITACT=1 and ILAM=1 . Default 0.0159 

Input Block 3: Integer Parameters. 

This block of input is read by the main program via a namelist called 

IINP . 

IMAX Number of vertical grid lines in the horizontal direction on the 
first grid. 1=1 corresponds to upstream infinity and I=IMAX 
corresponds to downstream infinity. For each grid refinement IMAX Is 
increased such that the new IMAX is two times the old value minus 
one. The limit on IMAX is 99. Default 13 

JMAX Number of horizontal grid lines in the vertical direction on the 

first grid. J=1 corresponds to infinity below the airfoil and J=JMAX 
is infinity above the airfoil. The same formula and limit that apply 
to IMAX also apply to JMAX. Default 7 

IKASE An integer number describing the case being computed. It is limited 
to a maximum of six digits and is printed at the beginning of the 
pressure printer plot for each grid. Default 100 

INV Parameter determining the program mode. It should be zero for 

analysis cases and one for inverse design cases. Default 0 

MITER Maximum number of interactions (complete relaxation cycles) allowed 
on the first grid. MITER is halved for each grid refinement. 

However, on the fourth grid, MITER is reset to 400. Default 1600 

NHALF Number of grid refinements. Default 2 

ITACT Viscous interaction control parameter. It should be set to zero for 
analysis cases without interaction and for design cases. It should 
be one for analysis cases with interaction. Default 0 
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ISKP2 Airfoil update control parameter for grid two. It should be zero if 
an airfoil shape update is desired on grid two every ten iterations. 
It should be one if an update is not desired until grid two solution 
is completed. Only used in the inverse design mode (INV-1). Default 
0 

ISKP3 Same as ISKP2 but for grid 3 

ISKP4 Same as ISKP2 but for grid 4 

ITERP Interpolation parameter for the design pressure distribution on grid 
four. If in the design mode the input pressure distribution for grid 
four is to be read as input data, INTERP should be zero. If it is 
desired to linearly interpolate the pressure distribution of grid 
three, it should be one. Default 0 

IREAD Starting solution control parameter. If IREAD is zero, the initial 
perturbation solution is assumed to be zero everywhere. It it is 
one, an initial solution is read as data. The latter would only 
normally occur when a user wished to restart a solution which had 
previously been saved. Default 0 

LP Relaxation cycle interval at which boundary layer details are 

printed. For diagnostic purposes suggest 50 or 100. For normal 
information purposes, suggest a value of 200. Default 1000 (no 
printout) 

IMASS Massive separation parameter. It should be one if the massive 

separation option is desired in analysis cases and is active only if 
ITACT is one. In inverse design cases (INV=1) it should be zero. 
Default 1 

I LAM Boundary layer parameter. If zero, boundary layer is computed as if 

all turbulent with transition at XIBDLY and XLBDLY. If one, boundary 
layer Is considered laminar- turbulent with natural transition. 

Default 1 

IPRT1 Print parameter. If one, perturbation potential values printed at 
the completion of each grid. Default 0 
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IPRT2 Print parameter. If one, x and y velocities at each grid point 
printed at the completion of each grid. Default 0 

KSEP Separated pressure distribution parameter. If zero, pressure in 

separated zone is assumed to be constant. If one, pressure in 
separated zone is considered variable. Default 1 

Input Block 4: Starting Solution (Optional). 

This block of data is read by subroutine VALUE only if the integer 

parameter IREAD has the value of one. 

P(I,J) Nondimensional perturbation potential at point I,J. Read by rows 

starting at J=JMAX down to J=l. Each row runs from 1=1 to I=IMAX and 
starts on a new line. Format 5E15.7 

PB(I) Nondimensional perturbation potential at point I on the y=0 (i.e. 

J=JB) grid line that is associated with the lower surface of the 
airfoil. Read from 1=1 to I=IMAX, Format 5E15.7 

Input Blo ck 5: Direc t Inverse Parameters (Inverse Design Mode only) 

This single line of input is read in subroutine COORD only when the 

inverse design mode is active (INV=1) . 

XI, X2 Same definition as in Block 3. However, when the inverse design mode 
is active, these values are read prior to the solution of each grid. 
This block corresponds to the first grid; and, thus, should always 
use X1=0 . 5 and X2=10000.0, Format 2F10.5 

Input Block 6; Airfoil Description 

This block of data is read by subroutine FOIL and describes the airfoil 

used in the analysis mode or the starting airfoil for the inverse design mode. 

NI The number of coordinate pairs used to describe the upper surface of 

the airfoil. Maximum value limited to 99. Format 15 
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XI (I) ,YI(I) Coordinate pairs describing the upper surface of the airfoil. The 
leading edge corresponds to XI=0.0 and the trailing edge is XI=1.0. 
The vertical ordinate, YI , is nondimensionalized by chord. Read 
starting with 1=1 to I=NI. Format 8F10.4 


DERIX , DERIY , DERFX , DERFY Parameters describing the leading and trailing edge 
of the airfoil. DERIX is dx/ds of the airfoil upper surface at the 
leading edge. It is usually zero. DERIY is dy/ds of the airfoil 
upper surface at the leading edge and it is usually 1.0. DERFX Is 
dx/ds 3 of the airfoil upper surface at the trailing edge. It is 
usually sufficiently accurate to use 0.0. DERFY is d 3 y/ds 3 of the 
airfoil upper surface at the trailing edge. It is usually 
sufficiently accurate to use 0.0. Format 4F10.4 


NIB The number of coordinate pairs used to describe the lower surface of 

the airfoil. Maximum value limited to 99. Format 15 


XIB(I) ,YIB(I) Coordinate pairs describing the lower surface of the airfoil. 
The leading edge corresponds to XI=0.0 and the trailing edge is 
XI=1 . 0 . The vertical ordinate, YIP, is nondimensionalized by the 
chord. Read starting with 1=1 to I=NIB. Format 8F10.4 

DERIXB , DERIYB , DERFXB , DERFYB Parameters describing the leading and trailing 

edge of the airfoil. DERIXB is dx/ds of the airfoil lower surface at 
the leading edge. It is usually zero. DERIYB is dy/ds of the 
airfoil lower surface at the leading edge and it is usually -1.0. 

e . 
of 

the airfoil lower surface at the trailing edge. It is usually 
sufficiently accurate to use 0.0. Format 4F10.4 

Input Block 7: Starting Airfoil Description (Optional) 

This block of data is read from subroutine FOIL and is only read if the 
integer input parameter IREAD is one. It effectively overwrites the 
information from input block six. 

YU(I) ,YL(I) , SLU(I) , SLL(I) Values describing the airfoil on the starting grid. 
YU(I) and YL(I) are the upper and lower surface ordinates, 
nondimensionalized by chord, at chord location X(I). SLU(I) and 


DERFXB is d 3 x/ds 3 of the airfoil lower surface at the trailing ed^ 
It is usually sufficiently accurate to use 0.0. DERFYB is d 3 y/ds" 
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SLL(I) are the upper and lower surface slopes at chord location X(I) . 
The X(I) values depend upon the size and spacing associated with the 
starting grid. The group of four values is read starting at the I 
value corresponding to the first point downstream of the leading edge 
(I=ILE) and ending with the point just upstream of the trailing edge 
(I=ITE) . Format 5E15.7 

DUPOLD ( I ) , DLWOLD ( 1J Values describing^ the boundary layer displacement 

thickness on the starting grid, DUPOLD(I) and DLWOLD(I) are the 
upper and lower surface displacement thicknesses corresponding to the 
chord location X(I) . These are read starting at the I value 
corresponding to the first point downstream of the leading edge 
(I=ILE) and ending with the point just upstream of the trailing edge 
(I=ITE) . Format 5E15.7 

Input Block 8 : Design Pressure Distribution (Inverse Design Mode only) 

This block of data consists of four sections which are only included In 

the inverse design mode (INV=1) . In that mode only the last three sections 

would usually be included. 

Section 1.- Starting solution design pressure distribution read by subroutine 
FOIL. This section would only be included if a design solution were being 
restarted (i.e. INV=1, IREAD=1, and MHALF=1) and it would only affect the 
first grid considered. 

CPU(I) Upper surface inverse region pressure coefficient values for the 

design case starting with 1=11, which is the first grid point after 
XI and ending with I=ITE, the grid point just upstream of the 
trailing edge. Format 8E10.3 

CPL(I) Lower surface inverse region pressure coefficient values for the 

design case starting with 1=11, which is the first grid point after 
XI and ending with I=ITE, the grid point just upstream of the 
trailing edge. Format 8E10.3 

Section 2.- This section reads in the starting and ending points of the 
inverse region from subroutine COORD and the inverse design pressure 
distribution from subroutine FOIL for the second grid. Used only in the 
inverse design case (INV=1) . 
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XI, X2 XI is the location where the direct calculation stops and the inverse 
calculation begins. Typically, it is slightly less than the third 
point from the leading edge or larger. X2 is the location where the 
inverse calculation stops. It should always be set to 0.5. Format 
2F10.5 

CPU(I) Upper surface inverse region pressure coefficient values for the 

design case starting with I— II, which is the first grid point after 
XI and ending with I=ITE, the grid point just upstream of the 
trailing edge. Format 8E10.3 

CPL(I) Lower surface inverse region pressure coefficient values for the 

design case starting with 1=11, which is the first grid point after 
XI and ending with I-ITE, the grid point just upstream of the 
trailing edge. Format 8E10.3 

Section 3.- This section reads in the starting and ending points of the 
inverse region from subroutine COORD and the inverse design pressure 
distribution from subroutine FOIL for the third grid. Used only in the 
inverse design case (INV=1) . The input variables and descriptions are the 
same as Section 2 above. 

Section 4.- This section reads in the starting and ending points of the 
inverse region from subroutine COORD and the inverse design pressure 
distribution from subroutine FOIL for the fourth grid. Used only in the 
inverse design case (INV=1) when ITERP is zero. Note that In references 2 
and 4, the use of grid four for inverse design is not recommended. The 
input variables and descriptions are the same as Section 2 above. 
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OUTPUT DESCRIPTION 

The printed output when the program is operated in the inverse design 
mode is identical to that described in reference 4. When the program is 
operated in the analysis mode with the massive separation option, the output 
for grids two, three, and four has the form shown below. Since the first grid 
assumes inviscid flow only, its printout only includes those portions 
associated with an inviscid solution. 

1. Heading (user supplied) 

2 . Case Number 

3. Mach number and angle of attack 

4. Case type callouts, i.e., 

INVISCID ANALYSIS CASE 

WITH LAMINAR TURBULENT VISCOUS INTERACTION 
AND MASSIVE SEPARATION 

AND VARIABLE PRESSURE IN SEPARATED REGION 

5. Input data in namelists FINP and IINP 

6. Coordinate System for current grid printed as I,X(I) followed by J, Y(J) 

7. Ordinates of the current airfoil displacement surface 

X- -Horizontal ordinate, where -0.5 is the leading edge and 0.5 is the 
trailing edge. 

YU- -Upper displacement surface ordinate 
YL- -Lower displacement surface ordinate 
UPPER SLOPE- -Slope of upper displacement surface 
LOWER SLOPE- -Slope of lower displacement surface 

8. Iteration history at ten-cycle intervals 
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CIR- -Circulation 


DPM- -Maximum <j> correction (absolute value) in the last relaxation cycle 
with the corresponding (I,J) grid location. 

NSSP- -Number of supersonic grid points 

DELTAY OR YUTE--In the design case, the change in YU(ITE) , the upper 
surface trailing edge ordinate, since the last surface update. 

Should go to zero if converging. In the viscous analysis case, it is 
the current value of YU(ITE); and it should approach a constant value 
if the solution is converging. Only changes after the first fifty 
cycles on each grid. 

SEP AT X--The current x ordinate value for the upper surface separation 
point (x — -0.5 is the leading edge and x = 0.5 is the trailing 
edge) . 

9 . Boundary Layer Information 

Every LP cycles results of the current boundary layer solution are 
printed first for the lower surface of the airfoil and then for the upper 
surface. In each case, the laminar solution (if ILAM=1) is printed first 
followed by the turbulent solution. 

9a. Laminar Boundary Layer Information (Printed every LP cycles) 

X- -Horizontal ordinate 
MACH #- -Local Mach Number 


CF- -Skin friction coefficient (the 0.1E11 initial value is arbitrary 
and should be ignored) 

D-STAR- -5*/ c > non-dimensional boundary layer displacement thickness 
D-THETA- -8/c , non-dimensional momentum thickness 
H -~S*/8, shape factor 
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RE -THETA- -Local Reynolds number based on B 


RE- STAR- -Local Reynolds number based on 5* 

TM- - 6 (du/ dx)/V , Pressure gradient parameter 

9b. Possible Laminar Boundary Layer Messages 

1. SEPARATION OCCURRED AT X — O.xxx -- gives x location where 
laminar separation occurred. 

. 2. SHORT BUBBLE FORMED? TRANSITION TO TURBULENT FLOW ASSUMED, 

X= 0 . xxx 

3. LONG BUBBLE? LAMINAR STALL MAY OCCUR, X - O.xxx BOUNDARY LAYER 
CALCULATION WILL BE CONTINUED AS TURBULENT BUT ACCURACY OF 
RESULTS IS QUESTIONABLE 

4. BOUNDARY LAYER CALCULATION COMPLETED? NEITHER SEPARATION FOR 
TRANSITION WAS DETECTED 

9c. Turbulent Boundary Layer Information (Printed every LP cycles) 

X- -Horizontal ordinate (-0.5 is leading edge, 0.5 Is trailing edge) 

i 

MACH #- -Local Mach number 

DLSTR--6 /c, non-dimensional boundary layer displacement thickness 
DEL - -6 /c, non-dimensional boundary layer thicknesses 
CUE--U e> transformed boundary layer edge velocity, U = (a <jo /a e )u 
USTAR- -Law of the wall parameter, SQRT(r w /p) 

USTAR**2- -Skin friction parameter, (r^/p) 

THETA- -0/c, non-dimensional momentum thickness 
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CDP- -Profile drag coefficient using Squire-Young approach modified 
according to reference 9 

CDW--Wave drag coefficient, not computed in this version of code and 
thus always zero 

CDTOT- -Total drag coefficient 

9d. Possible Turbulent Boundary Layer Messages 

1. SEPARATED CP IS O.xxx -- In the constant separated pressure level 
option (IMASS=1 , ITACT=1 ,KSEP=0) the pressure in the separated 
region is printed after each boundary layer calculation. If 
solution is converging, it should approach a constant value. 

2. USTAR2 (or USTCK) LT ZERO -- Indicates where in computation 
separation was first detected. See program listing for details. 

10. Final Boundary Layer Results 

YUORIG -- Original airfoil upper surface ordinate 
DU - - Smoothed upper surface displacement thickness 
SLU -- Slope of upper displacement surface 
YLORIG -- Original airfoil lower surface ordinates 
DL - - Smoothed lower surface displacement thickness 
SLL -- Slope of lower displacement surface 

11. Pressure Distribution on Airfoil 

CPU -- Upper surface pressure coefficient 
CPL -- Lower surface pressure coefficient 

12. Final Displacement Surface Information 
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YU -- Ordinate of upper displacement surface 


YL -- Ordinate of lower displacement surface 
SLU -- Slope of upper displacement surface 
SLL -- Slope of lower displacement surface 

13. Mach Chart 

The Mach numbers at the I,J coordinate points are multiplied by 100 and 
printed out in block form. The grid points "inside" the upper and lower 
displacement surfaces are indicated by zeros. Velocities (U,V) at the 
flowfield grid points may also be printed out using option IPRT2=1. 

14. Miscellaneous Information 

The normal force coefficient, CN, and the drag coefficient, WAVE CD, 
obtained by integration of the pressure distribution are printed out. The 
latter should theoretically be zero for subcritical unseparated cases, but it 
usually is non- zero due to mesh size and grid placement. Thus values of WAVE 
CD should only be used for comparison purposes. 

15. Printer Plot of Results 

U — Upper surface pressure coefficient 

L -- Lower surface pressure coefficient 
T - - Upper displacement surface 
B -- Lower displacement surface 

CPSTAR -- Pressure coefficient for local Mach number of one 

CLCIR -- Lift coefficient from computed circulation 

CL -- Lift coefficient from integration of pressure distribution 

CD - - Wave drag plus profile drag coefficient (accuracy depends upon 
value of CDWAVE) 

CMLE -- Moment coefficient about the leading edge 
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CDF -- Profile drag coefficient using Squire -Young approach modified 
for separation and compressibility 

CMC4 -- Moment coefficient about quarter chord point 

16. Miscellaneous Messages 

PROD .LE. ZERO 

This message indicates that the turbulent boundary calculation was 
unable to obtain an appropriate starting solution. As a result the 
flowfield calculations were continued with the displacement surface 
ordinates and slopes frozen at their last updated values. This 
situation usually occurs when transition takes place at the shock 
wave and the local Mach number immediately upstream of the shock wave 
is greater than 1.35. It probably indicates that at a minimum there 
is local separation in the vicinity of the shock wave. Since this 
phenomenon is not modeled in the present code the results obtained in 
such cases should be used carefully. Sometimes this problem can be 
"avoided" by computing an all turbulent case (ILAM=0) with transition 
at or immediately upstream of the shock wave. It should be noted 
that when this message appears , the subsequent computed values of CDF 
are in error. 


AA .LT. ZERO 

This message occurs when the local speed of sound is computed to be 
negative. It indicates that either the solution has become unstable 
or else the rotated scheme tried to use a point outside of the 
solution domain due to a supersonic point on the JMAX-1 row. Usual 
solution is to increase the supersonic damping, EPSS, and/or increase 
the stretching and/or size of the Y grid. When encountered, the 
computation is terminated and the current perturbation flowfield 
solution is printed for diagnostic purposes. 

17. Additional Output 

Note that additional output can be obtained for either diagnostic or 
analysis purposes by removing the C's from various commented print statements 
in the program. 
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TYPICAL RESULTS 


In the development of the present method, results have been obtained for 

NACA 4412 and NACA 0012 airfoils for freestream Mach numbers up to 0.6, angles 

of attack up to 18.5 degrees, and Reynolds numbers between three and nine 

million. These conditions were selected not because of limitations in the 

method but due to the availability of excellent experimental pressure 

22-23 

distribution data in those regions . Thus, the results presented here are 

only meant to be representative. 

Figure 4 and 5 compare results obtained with the present method with the 
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low speed experimental data of Pinkerton for a NACA 4412 airfoil at 6.3 
million Reynolds number. In figure 4 the experimental data have been plotted 
using the angle -of -attack correction suggested In reference 22. As can be 
seen, the theory predicts slightly larger lift coefficients than the data at 
the lower angles -of -attack . Whether this difference is due to an 
underestimation of the angle of attack correction, as suggested in reference 
24, or a problem in the theoretical model is unknown. In any event, the 
theory and the data are in excellent agreement between ten and fifteen 
degrees; and the present model reasonably predicts the location of maximum 
lift at an angle of attack of about 16 to 17 degrees. The theoretical model 
does, however, overpredict slightly (1.8 vs. 1.7) the maximum lift coefficient 
predicted by this forty- seven year old data. 

Figure 5 compares pressure distribution results obtained with the present 

method with the experimental data of Pinkerton at an angle of attack slightly 

below that corresponding to maximum lift. In this case, the corrected angle 

of attack was used in the computation, and the upper and lower surface 

boundary layers were assumed to be initially laminar followed by natural 

transition to turbulent flow. For this high lift case, the theory predicts 

that the lower surface remains entirely laminar and that the upper surface 

transitions at one percent chord followed by separation at 74.9 percent chord. 

As can be seen on the figure, the predicted pressure distribution is in 

excellent agreement with the data; and the pressure coefficient in the 

separated zone is slightly negative and constant. Experience indicates that 

for low speed cases better results are usually obtained using the constant 

pressure option (KSEP=0) for the separated zone. For this case, the 

theoretical lift coefficient was 1.69 while the experimental value was 1.68. 

The predicted profile drag coefficient was 0.0200, which is in reasonable 
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agreement with available measurements 

As Indicated previously, it is important for this type of method to 
include the effects of a laminar boundary layer. Figure 6 shows lower surface 


30 



laminar- turbulent transition point locations predicted by the present method 
for a NACA 0012 airfoil at Mach 0.3 and six million Reynolds number. As can 
be seen, for these flight conditions the lower surface boundary layer is 
predominantly laminar at low angles of attack; and at angles of attack above 
ten degrees it is essentially all laminar. Obviously, a method which only 
includes a turbulent boundary layer calculation and/or assumes transition near 
the leading edge might yield incorrect results for these flight conditions. 

Figure 6 also shows for the same conditions the predicted upper surface 
separation points determined by the TAMSEP method. Notice that no significant 
upper surface separation is predicted until about 12 degrees angle-of- attack. 
After that, as the angle of attack increases, the beginning of separation 
moves forward on the upper surface until more than half of the airfoil 
experiences separated flow at about 16 degrees. For this case, maximum lift 
occurs around 14 degrees . 

The experimental data from the Low Turbulence Pressure Tunnel which are 
presented in figures 7, 8, and 9 were obtained from Mr. C. L. Ladson of the 
NASA Langley Research Center in Hampton, Virginia. Figure 7 compares the 
predicted lift coefficient variation with angle of attack for a NACA 0012 
airfoil at the same nominal conditions (i.e. Mach 0.3 and six million Reynolds 
number) with experimental data obtained in the Low Turbulence Pressure Tunnel 
at NASA Langley. These data were obtained on clean airfoils without the use 
of trip strips, and thus the theoretical results were obtained using the 
laminar natural- transition turbulent model (ILAM=1). At the lower angles of 
attack the agreement is quite good. However, at the onset of trailing edge 
separation, around twelve degrees, the theoretical lift curve exhibits a 
"kink" accompanied by a slight decrease in lift coefficient. This "kink" is 
often observed in the theoretical lift curve when separation is first 
predicted and is due to the fact that separation is usually first detected on 
the coarse grid. On that grid, the first point forward of the trailing edge 
is about 93% chord; and, thus, when separation is first predicted on the 
coarse grid the separation point moves forward from the trailing edge to at 
least the 93 percent point, with the result that the amount of separation is 
overpredicted and the lift at that condition is underpredicted. Since the 
model only permits the separation point to move forward, this effect is 
maintained throughout all the grids at that angle of attack. However, as the 
angle of attack is increased this effect disappears. Thus, at low and medium 
freestream Mach numbers, the lift coefficients predicted at angles of attack 
just above the onset of trailing edge separation are usually slightly low. 

As can be seen on figure 7, the lift coefficients predicted at the higher 
angles of attack are, at least for this case, in reasonable agreement with the 
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experimental data. Notice that for this case the experimental data indicate 
an apparent maximum lift around 14 degrees followed by a decrease and then an 
increase. The theoretical lift based upon the calculated circulation predicts 
a maximum lift coefficient of about 1.40 at 14 to 16 degrees. On the other 
hand, the theoretical lift based upon integration of the pressure 
distributions indicates a maximum lift at 14 degrees, which nominally agrees 
with experimental data. This slight divergence between the values predicted 
by circulation and those by pressure integration is, based upon experience, an 
excellent indication of the maximum lift location at low and medium speeds. 

The reason for this statement will be evident when the pressure distributions 
for these cases are discussed. 

The predicted variation of drag with lift for this case is compared to 
experimental data at the higher lift coefficients on figure 8. Again the zig- 
zag in the theoretical curve corresponds to the similar phenomena on the lift 
versus angles-of-attack plot and is due to the overprediction of the size of 
the initial separated zone. Nevertheless, the agreement between the 
theoretical predictions and experimental values, particularly near maximum 
lift, is good and should be acceptable for applied engineering calculations . 

Figures 9(a-c) compare presure distributions obtained with the present 
method with data obtained by Ladson in the Low Turbulence Pressure Tunnel at 
three different angles of attack. The first corresponds to an unseparated 
flow situation, the second is near maximum lift, and the third is for an angle 
of attack above the maximum lift condition. Since evidence at medium 
freestream Mach numbers and higher indicates that the pressure variation in a 
separated zone is not constant, these and subsequent cases were all run using 
the variable pressure option (KSEP=1) . As can be seen on figure 9(a) the 
theoretical pressure distribution for the unseparated case is in excellent 
agreement with experimental data. For this case, the theoretical lift and 
drag coefficients were 1.27 and 0.0116 while the corresponding clean airfoil 
experimental values were 1.23 and 0.0123. 

For the case near maximum lift, figure 9(b), it should be noted that some 
supersonic flow exists over the upper surface of the airfoil in a very small 
region near the leading edge since the critical C p for this case is -6.89. In 
addition, the theoretical method predicts upper surface separation at 74.9 
percent chord and boundary layer instability on the lower surface at about 80 
percent chord. However, due to the favorable pressure gradient, the lower 
surface boundary layer never transitions. For this case the theoretical lift 
coefficient of 1.39 coincides with the experimentally measured value, and the 
two pressure distributions exhibit reasonable agreement. 
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At an angle of attack greater than that corresponding to maximum lift, 
the flow about an airfoil is typically characterized by a large region of 
unsteady separated flow; and a steady state solution method such as the 
present one is not really applicable. Thus, the apparent lack of agreement 
between the present steady theory and the experimental measurements shown on 
figure 9(c) is not surprising. Nevertheless, the general pattern of the 
pressure distribution including the existence of a large separation zone is 
predicted; and the predicted lift coefficient of 1.44 based upon circulation 
is in surprising agreement with the experimentally measured values of 1.437. 
However, careful examination of the solution indicates that it is not 
completely converged, that the theoretical lift may be oscillating slightly, 
and that the lift based upon pressure integration computed at the end of a run 
using default parameters is only 1.27. Interestingly, results using the thin 
layer Navier Stokes equations 20,21 for a similar case (NACA 0012, Mach 
No. =0.3, Reynolds = 1 million, angle of attack = 18 degrees) indicate a lift 
coefficient varying with time from 0.65 to 1.6 with a Strouhal number of 0.1; 
and both the theoretical and experimental pressure profiles shown on figure 
9(c) are representative of those computed at various times in the cycle with 
the thin layer Navier Stokes model. Thus, the present theoretical result is 
representative of the type of pressure distribution and lift which might exist 
for this condition. 

Another interesting feature associated with the results shown on figure 
9(c) is that the drag coefficient predicted using the method of reference 6 
(i.e. CDF) was 0.0604 while that predicted using the method of reference 9 
(i.e. CDP) was 0.1190. Normally, these two values are in good agreement with 
each other. Apparently, in the present method when the maximum lift condition 
is exceeded the solution becomes oscillatory and not completely converged, the 
lift computed by pressure integration diverges from and is lower than that 
from circulation, and the CDF and CDP values differ significantly. It is 
believed that these three items can be used to determine for medium Mach 
numbers the angle of attack corresponding to maximum lift. 

Figure 10 shows for the same NACA 0012, Mach 0.3, Reynolds number 6 
million case, the displacement surfaces predicted by the present method for 
the upper surface region between 70 percent chord and the trailing edge at 
various angles of attack. Below 11.13 degrees, where the flow is unseparated, 
the displacement thicknesses in the trailing edge zone are relatively small. 
However, with the onset of separation at 12.09 degrees the thicknesses begin 
to increase rapidly, and the displacement surfaces take on shapes 
characteristic of the flow over a stalled airfoil. It should be noticed that 
at an angle of attack of 16.12 degrees, the displacement surface starts to 
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curve back towards the freestream angle prior to the trailing edge. It Is 
believed that the displacement surfaces shown on this figure have the correct 
behavior and are an adequate engineering representation of the real flow. 

Obviously, it would be desirable for the present method to accurately 
model transonic flows with and without significant trailing edge separation. 
Consequently, predictions obtained with TAMSEP have been compared with data 
obtained in the NASA 8 Ft. Transonic Pressure Tunnel by Harris 22 . These data 
probably represent the best high lift transonic experimental airfoil data 
available today, and it has been used by many investigators. However, in most 
cases previous studies have compared results using the tunnel geometric angle- 
of-attack values and have ignored the known corrections associated with these 
data. Since such comparisons could lead to erroneous conclusions, the present 
results were obtained by matching the tunnel Mach and Reynolds numbers and 
using the corrected angles of attack suggested by Harris. 

Figure 11 shows a transonic separated flow result compared with data 
obtained by Harris 22 . While the pressure distribution shown was obtained 
using the laminar- turbulent boundary layer option, indistinguishable results 
were obtained assuming transition at 6 percent chord. In the actual 
experiment, boundary layer trip strips were located at 5 percent chord. For 
this case, the present method predicts upper surface separation at 87 percent 
chord and lower surface transition very near the trailing edge. While the 
experimental shock location is slightly forward of the theoretical value and 
while the experimental pressures in the trailing edge region are slightly 
lower than those predicted by the theory, the overall agreement is probably 
acceptable for engineering studies. For this case, the measured normal force 
coefficient was 0.994 while the predicted lift coefficient was 0.981. 

A lift versus angle-of-attack curve typical of those predicted by the 
present method is compared with experimental data on figure 12. For this 
freestream Mach number of 0.5, significant transonic flow accompanied by a 
strong shock wave is present on the upper surface at angles of attack greater 
than 6 degrees. As can be seen, the theoretical prediction, which was 
obtained assuming transition at 6 percent chord, agrees well with the 
experimental data up to about 7.35 degrees. Above that angle of attack, the 
present method predicts a maximum lift coefficient of 1.09 at 8.24 degrees 
with no trailing edge separation. At 9.21 degrees the theory predicts a 
decrease in lift coefficient to 0.982 with upper surface separation at 87 
percent chord. The experimental data, however, indicates that maximum lift is 
1.02 at 9.21 degrees and that trailing edge separation probably started at 
about 7.5 degrees. Examination of the theoretical results at 7.35 degrees 
reveals that the local Mach number immediately upstream of the upper surface 
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shock wave is 1.42. Such a strong shock wave, whose strength increases with 
angle of attack, should induce significant shock boundary layer interaction 
which, unfortunately, is not modeled in the present theory and code. 

Therefore, in this case the differences between the theoretical results and 
the experimental data near maximum lift are probably due to shock boundary 
layer interaction and its subsequent effect on boundary layer growth and 
trailing edge separation. Nevertheless, the present model and code does give 
a reasonable indication of the location and magnitude of the maximum lift 
coefficient. 

A similar lift curve for a freestream Mach number of 0.55 is shown on 
figure 13. As before, the theoretical results were obtained using an all 
turbulent boundary layer; and in this case no upper surface trailing edge 
separation was detected until an angle-of -attack of 8.266 degrees. The 
maximum lift coefficient was computed to be 1.02 at 7.29 degrees as compared 
to the experimental values of 0.983 and 8.266 degrees. For this case 
significant transonic flow existed at all angles -of- attack above 4.5 degrees 
and by 7.29 degrees the Mach number at the upper surface shock wave had 
increased to 1.50. At 9.33 degrees a complete solution could not be obtained 
with the present method. On the medium grid the upper surface flow separated 
at 87 percent chord, and on the fine grid the separation point moved forward 
to the shock wave at about 18 percent chord and the solution failed. Quite 
obviously significant shock boundary layer interaction exists at these high 
angles of attack and the decrease in lift or stall is probably due more to 
shock induced separation than to the onset of significant trailing edge 
separation. 

Nevertheless, the present method can be used to estimate reasonably 
accurately the occurrence of this situation. As can be seen on figure 13, the 
method predicts reasonably well the magnitude of the maximum lift coefficient 
and is conservative as to the corresponding angle-of -attack location. In 
addition, by noting the mechanism of code "failure”, in this case sudden 
separation at the shock wave, a user can probably determine the type of stall 
phenomena present. 

Theoretical pressure distribution results are compared with experimental 
data for a freestream Mach number 0.6 case in figure 14. This case at a 
corrected angle of attack of 5.59 degrees and three million Reynolds number Is 
significant for a variety of reasons. First, it is an example of a transonic 
case with a strong upper surface shock wave. Second, the flow is unseparated 
and the data should serve as a good test of the present method for a situation 
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without separation. Finally, this case has also been solved by Anderson 

20 

et al. using both an Euler boundary layer method and a thin layer Navier 
Stokes method. 

It should be noticed that the present results, like those of Anderson, 
agree very well with the experimental pressure distribution with respect to 
shock location and pressure levels. Also, for this case the experimental lift 
coefficient was 0,781. The present TAMSEP method predicted 0.809, the Euler 
boundary layer method of Anderson yielded 0.804, while his thin layer Navier 
Stokes result was 0.793. Obviously, the present method Is capable of yielding 
excellent results that are in agreement with experimental data and other 
analytical methods. However, since it is a non- conservative full potential 
method, it should obtain such results with accurate shock wave locations more 
easily and faster than other more complicated methods. 

Another lift versus angle-of-attack comparison is shown on figure 15 for 
the NACA 0012 at a freestream Mach number of 0.6 and a Reynolds number of nine 
million. As can be seen the highest theoretical point plotted is at a lift 
coefficient value of 0.945 and an angle of attack of 6,392 degrees. At 7.348 
degrees the boundary layer solution became frozen (i.e. PROD. LE. ZERO) on the 
medium grid due to shock boundary layer interaction, and on the fine grid the 
flow separated at the shock wave and a converged solution was not obtained. At 
8.371 degrees the flow separated at the upper surface shock wave on the coarse 
grid, and again a converged solution was not obtained. For these cases, 
computed local Mach numbers in the vicinity of the shock wave were as high as 
1.56. Thus, the theory indicates that above 6.392 degrees significant shock 
boundary layer interaction probably accompanied by separation exists and that 
the maximum lift coefficient occurs at that angle-of-attack. As can be seen 
on the figure, the magnitude of the predicted maximum lift coefficient is in 
good agreement with the experimental data; although the angle-of-attack 
location is again conservative. 

Based upon the results presented in this section, it is believed that the 
present method and code can be used at low and medium Mach numbers to 
accurately predict lift and pressure distributions at angles of attack up to 
that associated with maximum lift. At transonic speeds, the method should 
give good results for unseparated flows and for flows having trailing edge 
separation without significant shock boundary layer interaction. Thus, at 
transonic conditions, the method is probably currently limited, for accurate 
results, to Reynolds numbers of three million and higher and to local upper 
surface Mach numbers less than 1.4 to 1.45. In addition, it should yield 
reasonable estimates for the maximum lift coefficient at transonic speeds 
while being conservative as to the corresponding value of angle of attack; and 
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the method should indicate the onset of significant shock boundary layer 
interaction. 


CONCLUSIONS 

A direct- inverse technique based upon a nonconservative full potential 
inviscid method, a Thwaites laminar boundary layer technique, and the Barnwell 
turbulent momentum integral method has been developed. This method is 
suitable for predicting the subsonic and transonic flowfield about airfoils 
having trailing edge separated flow. Extensive comparisons with experimental 
data indicate that the method should be a useful tool for applied aerodynamic 
engineering analyses. In addition, it is believed that the range of 
applicability of the method could be extended significantly by the addition of 
a shock boundary layer interaction mode. 
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Direct solution with 
boundary layer interaction 



"| C p constant 
► boundary 
J condition 


Figure 1.- Problem formulation. 
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Figure 2.- Iteration and interaction procedure used in TAMSEP. 
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Figure 3.- Subroutine tree for TAMSEP. 
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Figure 4.- Comparison of predicted lift coefficient with experimental data 
NACA 4412, M - 0.15, Rn - 6.3 x 10 6 . 
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Figure 5.- Comparison of theoretical pressure distribution with experimental 

data. NACA 4412, M =0.15, Rn = 6 . 3 x 10 6 , a = 15 . 3 deg. 
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Figure 6.- Lower surface transition and 
NACA 0012, M - 0.30, Rn - 6 
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Figure 8.- Comparison of theoretical lift and drag with experimental data 
NACA 0012, = 0.30, Rn = 6 x 10 6 . 
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(a) a = 11.3 degrees. 

Figure 9.- Theoretical and experimental pressure distribution comparisons. 
NACA 0012, = 0.30, Rn^ - 6 x 10 6 . 
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(b) a = 14.3 degrees. 
Figure 9.- Continued. 
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(c) a - 16.1 degrees. 
Figure 9.- Concluded. 
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Figure 11.- Theoretical and experimental pressure distribution comparison. 
NACA 0012, - 0.5, Rn tt = 6 x 10 6 , a = 9.25. 
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Figure 12.- Comparison of predicted lift coefficient with experimental data. 
NACA 0012, M = 0.5, Rn - 9 x 10 6 . 
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Figure 13.- Comparison of predicted lift coefficient with experimental data 
NACA 0012, = 0.55, Rn^ = 9 x 10 6 . 
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Figure 14.- Theoretical and experimental pressure distribution comparison. 
NACA 0012, = 0.6, Rn a - 3 x 10 6 , a - 5.59 deg. 
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TEST CASE, NACA 0012, MACH 0.5, RN 6 MILLION, ALPHA 9.25 DEG 
&FINP M=0 . 50 1 , ALP=9 . 250 , RN=6 . 00+06 , 

EPSS= 1 . , XPC = . 5 , &END 
&IINP IKASE=509 ,NHALF=2 , SEND 
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PLANE-FREE STREAM FROM TOP 


1*2, 12 T OF TO SOT TOM 
J* 7 , 4 LEFT TO RIGHT 


SO 49 48 81 S3 

2 S IB 0 10 1 44 

43 43 0102 41 

4 S 41 0 77 17 

41 17 07114 

4 4 47 O BS I 1 

44 47 0 It BB 

48 48 0 87 BB 

44 4 S 0 42 62 

36 33 0 41 SO 


NORMAL FORCE COEFFICIENT CN BY INTEGRATION « 1.2181 

WAYE CD • 010318* 

THEORETICALLY ZERO FOR SUICR1TICAL CASES 
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CPSTAR • -2.1317 

« 0 103114 C ML E * 


CLC I R ■ 
•0.3013 CDF 


1 .34 11 

O . O 


TEST CASE, 


HACA 0013, MACH OS. RN S 
CASE NUMBER 
MACH NO IS 0 . SOIANCLE 


MILLION, ALPHA 9.2* DEC 
KOI 

OF ATTACK IS 9.210 DECREES 


IFINF 

M * .*00999917 ,W« 

1 . OOOOOOOO ,**• 

3.19999919 .RN* 

RDELFN* . 12*000000 

xliep* .*00000000 

BEND 

1 1 I NP 

I MA X * 2 » , JMA K 

O, 1SKR3* 

1 , IRRT 1 * 

BEND 


Irni'JISiiJrTuiiuiSIt VISCOUS INTERACTION 

InO VARl ABLE*PRIISURE IN SEPARATED RlCION f p IS , 

| 9BB9BBB 1 ,*»• *00000000 .**« 1 LSSJS*, • o/i M llilitMa ! « !’ 1 4999997 S ,»• 

1*00000 10 ,S4. 2 00000000 ' C 0 “ V ^ ° ° J ' C DC ORR • ” ."OEM 2*0000000 
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470000031 
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SOI . I NV « 
© , I TERR • 

0 . KlIFi 
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0 . I Rf AD* 


tOO.NHALF* 2 • ITACT* 

o , L P > 1 OOO , IMASS • 


l . ISKR2 
I.ILAMi 


X - Y GRID SYSTEM 


3 •© . 39729 *0 I 

• - 0 30*19*00 
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2 -0.111 11*00 

• 0. 45929-01 
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t • 0.1 24 7E*O0 

2 1 O . 4 R OOE *00 

3 .0 42* 1 **00 
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4 -O 14719*00 

1 O -0.1 1 S1E *00 

IE 0.1 1S1E *00 
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4 - O . 24*09+00 
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s - O . 4*001+00 
n -O 12 4 7 9* 00 

17 0.241*9*00 

23 0.14109*01 
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11 O . 42S 19*00 


1 - O 13079*00 

12 - 0 . I2AOE - O 1 
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• -0.1S92E-01 

12 0 . 9 1 B 19*00 


7 -0.37099*00 

13 0.0 

II 0 S709C*00 

7 *0 . 92 1 39 - 07 
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AIRFOIL COORDINATES 

* YU Yl UPPEJl SLOPE LOWER SLOPE 
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-O. 1**77 O 0*9*1 -00519* -0 005*8 0.005*5 
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NSSP t 

2 

OELTAY 

OR 

YUITE 

a 

O . 0392 

SEP 

AT 

X • 

0 37 1 

NS5P ■ 

2 

OELTAY 

OR 

VUITE 


O 0313 

SEP 

AT 

X • 

0 37 1 

NSSP * 

2 

OELTAY 

OR 

YUITE 


O 0393 

SEP 

A T 

X ■ 

0 3 7 1 

NSSP * 

2 

OELTAY 

OR 

VUITE 

a 

O 0314 

SEP 

AT 

X ■ 

O 37 1 

NSSP > 

2 

OELTAY 

OR 

YUITE 

• 

O 0305 

SEP 

AT 

X « 

0 3 7 1 

NSSP • 

2 

OELTAY 

OR 

YUITE 

« 

0 0315 

• CP 

AT 

X * 

0 37 1 

NSSP * 

2 

OELTAY 

OR 

VUITI 


0 0311 

sep 

AT 

X • 

0 37 1 

NSSP * 

2 

DE L TAY 

OR 

VUITE 

t 

0 0317 

SEP 

AT 

X ■ 

0 3 7 1 

NSSP * 

2 

DEL TAY 

OR 

YUITE 

a 

0 .0317 

ser 

AT 

X » 

0 37 1 

NSSP < 

2 

DEL TAY 

OR 

VUITE 

• 

0 0311 

SEP 

AT 

X * 

0.371 

NSSP * 

2 

OELTAY 

OR 

YU ITE 

• 

0.0311 

SEP 

AT 

X I 

0 37 1 

NSSP ■ 

2 

OELTAY 

OR 

YUITE 

a 

0 0311 

SEP 

AT 

X * 

0.371 
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OWGINAL PAGE IS 
OF POOR QUALITY 


BOUNDARY LAYER ANALYSIS FOR REYNOLDS NUMBER OF O.BOOE*07 


X 

0 . <»POO 
O 430S7 
0370** 
0 3 0 9*2* 
0 . 2*952 
O 1**77 
O 12470 
006240 

0 o 

00*240 
O 12*70 
O 19677 
0 24*52 
O 30*42 
0 37051 
0 4 3067 
O 41000 


YUOR I G 
0.01704 
O 04070 
O 05112 
0 05*77 
O 0**45 
O 05 • *■ 
005**5 
O 05 642 
005294 
0049*0 
004355 
0 037* 1 
00317* 
002*14 
001*0* 
0.0105* 
O 002*5 


DU 
O . O 

0.0001* 
0 . 0003* 
O . 00070 
0.001 lO 
0.00151 
0.001*4 
0.00240 
0 002* I 
00034* 
0 , 004 1 • 
0 005 1 0 

O 001*4 
0010*3 
0 01732 
0 02*45 
0 03*47 


SLU 

O . * 105* 
O . 19721 
0 13970 
0 . 0**47 
003171 
O 00017 
002318 
004025 
005441 
O .0*549 
O 07375 
0079*5 
00*235 
O 02754 
O 00131 
0037*3 
O 09 *03 


YL OR! C 
*0 01704 
-O . 04070 
•00*112 
' O 0* *77 
-O 06*45 
-O .05*** 
•0.05*16 
- O . 05642 
•O . 052*4 
• O . 04660 
•0.043*5 
• O . 037* 1 
•0.0317* 
-002514 
-0,01 *08 
•0.0105* 
- 0 002*5 


DL 
O . O 

O OOO 1 « 
O . OOO I * 
0 00022 
O 0002* 
O . 00030 
O 00035 
O 0003* 
O 0004 2 
O . O 00 4 * 
O 00050 
O 000*3 
O 0005* 
O 000*1 
O . 00083 
0.00109 
0.00135 


CP BY CENTRAL 
* CPU 

•0 490 

-0 431 
•0.371 
•6.310 

• 0 24* 

*0 1*7 

• O 12 5 
- O 0*2 

O . O 
O 0*2 
0 12* 
0.1*7 
O 249 
0 3 10 

0 37 1 

0 4 3 1 

0.4*0 
X 

-O 4*000 

• O 4 30*7 
•O 37058 

* O 309 9 2 
- O 24*52 
- O 1 • *77 
-O 12*70 
-00*240 

O O 

005240 
0 12470 

0 15177 

024*52 
0 309*2 
O 3705* 

0 43017 

0 49000 


D I PPCRCNCIS 

CP L 
1.049 
0 5*1 
0.350 
0 229 

0 160 
O 112 
0.0*1 
O 0*0 
0 047 

0 03* 
0 03 3 
O 030 
0 030 

0 032 
0.035 
0 04 3 
0 115 

V L 

-001703 
•O . 0405* 

- 0 05 1 30 

- 0 05 *• • 

*00*971 

• 0 . 0902 9 
•00*920 
•005991 

- 0 OS 339 

• 0 04 90* 

•0.04405 
-003*45 
•0 03234 

•0025*0 

- O 0 11*1 

• © 0 11*9 

-O 00402 


001703 
0.040** 
0.01111 
005748 
0.0*054 
0 0514* 
O . 0*075 
0 01117 
0.055*5 
005209 
0.04774 
0 04303 
0 Oil * 1 
0 03577 
0 034*1 
003*13 
003992 


- 2 722 
-2**7 

- 1 .71* 

- 1 .313 

- 1 075 

-0 **1 
•0.717 
-0.577 
*0.454 
-0.340 
' 0.220 
*© . 112 

O 002 
0 07* 
0.107 
0.157 
0 19* 


SLU 

0*1059 
0 19721 
O 13*70 
0.0*947 
003171 
O 000*7 
0 0221* 
0.04021 
00*444 
O 0**49 
O 0737* 
0.07111 
0.0*23* 
0.02754 
0.00131 
0037*3 
00*903 


5 L L 

0*1058 
O 1 9 5 6 2 
O 1 35 7* 
0.0833* 
0 02*79 
0 005 4 2 
002911 
0 04735 
OOI249 
0.07415 
0 0*551 
0094*3 
O 10300 
011012 
0 1 195 4 

O 12449 
O 13412 


MACH CHART IN COMPUTATIONAL 


SL L 

0.91058 
O 19652 
O 13575 
0.05325 
0.02579 
O . 00542 
0.02*56 
O 0*739 
00*249 
O . 0749* 
0.0*56* 
0.094*3 
O . 10300 
O 11012 
O . I I 154 
© 12449 

O 13412 


PLANE-FREE ^TRMM F^OM TOP 


I • 2 , 24 TOP TO BOTTOM 
j* 2,12 LIFT TO R I OMT 


*0 *0 90 60 50 


43 3* 31 24 14 
41 43 39 39 33 
49 43 40 39 39 
49 43 42 42 43 


4 7 4 7 47 49 49 

47 47 41 4* 49 

47 47 4* 4* 49 

*7 47 47 47 47 


10 60 60 50 50 90 

4* 49 60 50 51 51 


0101 97 7* 70 99 

0109 94 72 13 61 

O *• 79 71 13 9* 

O 91 79 9* 13 5* 

O 7* 72 *7 92 55 

0 71 *9 II 51 5* 

0 97 99 93 *0 9* 

O 14 *3 *1 59 9* 

0 |1 *0 69 56 5* 

O 5* 5* 57 If 64 

O 95 9* Cl 95 54 

O 13 53 54 54 64 

O 50 II 52 53 53 

0 41 60 It 12 53 

O 47 49 10 52 S3 

0474*605152 
O 4* ** 50 52 54 

49 49 49 50 50 *1 


NORMAL FORCE COEFFICIENT CN BY INTEGRATION 
WflVI CD • 0.057*24 

THEORETICALLY ZERO FOR SUBCNJTICAL CASES 


O *4*7 
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CASE NO 101 

-O MOO U TB L 

• o <300 U TIL 

* O . 3700 U TIL 


-O 2100 
•0.3300 
• 0 1700 

- O 1100 

- o osoo 
o O 1 oo 

0,0700 
O 1 300 
O 1 IOO 
0 2 5 00 

0.3100 


u 


T I l 

U TIL 

U TIL 

U TIL 

U T I 

U T I 

U T I 

U T ■ 


U TLB 
TL I 
TILU 


O 3900 

o tipo 

O 4100 


TILU 

TIL U 
TB L U 


• 2 117 


CL 


0 4* IS CD 


- 1 . 892 -1.301 

MISSURE COEFFICIENT 

CPSTA* « *2 1217 

0 079229 CMLt « 


•0 . 909 >0 & 17 -0.121 

CLC 1 I • 0 . 1203 

0.1997 CDF ( 0.010901 CHC 4 


0.299 0.994 


0.0999 


HIT CASE. 


NACA 0012, MAC M 0.1. Ill « 

CA9I MUMB1* 

MACH NO 19 0 90IANCL9 


MILLION, ALPHA 

BOB 

OP ATTACK IS 


9 29 DKG 
• 210 DKCftCIt 


AFINP 


1 OOOOOOOO , K4< 

3 99999949 ,NN« 

IDfLFM* .129000000 
XL SIP * . 900000000 

• INO 
1 I I NP 

IMAX- 4 9 , JMA X 

0 , I 9 K P 3 * 

1 , I BIT 1 * 

IIND 


IMVItCIO ANALYSIS CASK 

WITH L AM I MAI TUIBULINT VISCOUS INTIIACTION 
AND MASS 1 VI SIPAIAT ] ON 

AND VAIIASLC PllttUII IN SIPAIATIO IIC10N 


. 4900000 
9000000 . 
.91* 

,xpc* 


tit , X 1 • 

10 ,94* a 
00 .XIIDLV 
4000000991-02 
900000000 


391977944 .13* . 900000000 .ALP* .111442971 

OOOOOOOO , CDNV * 1 000000 341 - 09 , A 1 ■ 2499B9II2 

* *.439919919 .CII* 490149929 .COCONR* .O 

Ml*. 431999999 ,ICPI» 199919999 r CP9* 

, X L 10 L V * - 439999999 , RL*X< 1.00000000 . IADUS 


.BPS* O , IPS)* 

.42- . 149999979 .A3* 

. BOIL • . 290000000 

399999979 .XMON* 470000029 ~ 

* . 1990000091-0 1 


29 , I KASE • 
1 SKP4 • 

1 PRT2 » 


• 09 , INV * 
O , ITIIP * 

0 , KBIP* 


0 . MI Til • 
O . I II AD- 


400, NHALP* 2.ITACT- 1 , 1SKP2- 

O.LP* IOOO, 1MAJI- l.ILAM* 
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X - V GRID S VITIM 


2 - 0 . 10271+0 1 
• -O . 13171*00 

M * O . 34031 + 00 
20 - 0 . 1 1111*00 
24 031211-01 

32 0.2 1 7 7 ■ ♦ OO 

31 O. 40071*00 
44 0.11 111*00 

3 -O. 1 411I+01 
a - o ii 4 *i+oo 

14 O 32311-01 
20 0.32011*00 


3 -0.33721*01 

• -0.41001*00 
11 -0.30141*00 
21 -O 12*71*00 
27 O . 13401-0 1 
33 0.24111*00 
31 0.43071*00 
41 0.14101*01 

3 - 0 . 11 111*00 

• - 0 14201*00 
11 O . 11121-01 
21 0.42111*00 


4 -022711*01 

1 0 -O 41041*00 
If -O. 27121*00 

22 -0.13171-01 

21 0.13171-01 

34 0.27121*00 

40 O 41041*00 

41 0.22711+01 

4 -0.19311 +00 

10 - 010 111+00 
I I O 1 0 1 11 + 00 

23 O . 113*1+00 


1 -0.14 101*01 

11 -0 . 43071+00 

17 -0 24111+00 

23 -0. 12401-01 

21 0.1 2471+00 

31 0 . 30111+00 

*1 0 41001+00 

47 6 31721+01 

1 -0.42111+00 

11 -0.31121-01 

17 01 420E+00 

23 0111 11+00 


• - 0 . 11 111+00 
12 *0. 40071*00 

11 -0.31771*00 

24 -0.21311-01 

30 O 11131+00 

31 0.34031+00 

42 0 13171+00 

41 0 10271+01 

1 -0.32011+00 

12 -O. 32311-01 

II O . 11111 + 00 
24 0 . 11111 + 0 1 


7 -0 14711+00 
13 -0 370*1+00 
11 -0 . 11111+00 
21 0 0 

31 O 1 l 1 11+00 
37 0.37011+00 
43 0.34711+00 

7 -0 24101+00 
13 -0.13131-07 
11 0 . 24101+00 


•O 41000 

• O 4 104 4 

• O 43017 

• O . 40072 

-O 37011 
-O 34021 
•O . 30112 
-O . 27123 
- O . 24112 
-0.21711 
-O 1*177 
-© 11177 

- O 12470 
- O 01317 

• O 01240 
-0.03121 


vu 

0.01704 
003212 
004011 
0.04 120 
0.01111 
o . 0144* 
O . 05741 
0 . 0*10 1 
O 01014 
0.01101 
001141 
0 01114 
0 . 01071 
o out 1 
0 . 01112 
0.01734 
0.01111 


AIRFOIL COORDINATES 

VI UPFIA SLOPE 


-0.01704 
-O . 032 1 2 
-0.04014 
- 0 . 04104 
- O . OS I 30 
-0.01411 
•O . 01119 
•O . 01131 
-0 0*171 

- o oiooo 
■o.oioii 
•0 01171 
'0 01120 
- o oiioo 
• 0.01111 
- 0.01101 
•O 01331 


O . 1 ion 

0 31134 

O 11721 
O 11141 

O . 13170 
O . 10401 
0.01147 
O 01001 
003171 
0.01411 
O 00017 

- 0 0 1074 

-002211 
•0.03120 
-004021 
-o 04731 

- 0 0*44 1 


0.03121 

0 . 06397 

-0.01121 

-o mil 

0 

0 01240 

0 . 05309 

-0 . 04909 

-0.OI44I 

0 

0 . 01317 

O 04992 

-0 04111 

- 0 019 12 

0 

0 . 12470 

0 0*774 

-0 . 04401 

*0 07371 

0 

O . 11977 

0 0413* 

•0.04121 

-0.07120 

0 

O . 1 1177 

0.04303 

-003444 

-0.07111 

0 

021749 

0 0*012 

- 0 03*39 

- o o«i*o 

0 

0.24112 

0 . 0311 1 

-0 03234 

- 0 01235 

0 

027123 

0 037 I 1 

-0 02107 

-0.04494 

0 

0 30112 

0 . 03177 

-0 02*40 

- 0 027*4 

0 

0 . 34021 

0 03*37 

•0 02231 

-001311 

0 

O . 37011 

0 03419 

-0 . O i 11 1 

0.00131 

0 

0 . 40072 

0 03111 

•0.01130 

0 01917 

0 

0 43017 

003113 

-0.01111 

0 03713 

0 

O 4 104 4 

0.03103 

-O 00711 

0 OB343 

0 


O 4 B OOO 
I THAT I ON 
ITERATION 
ITERATION 

iteration 

I TERAT I ON 
I TERAT I ON 
1 TERAT ION 
I TERAT I ON 
ITERATION 
1 TERAT ION 


O . 03112 
10 C IR • 
20 C I R • 
30 C 1 R • 
40 C 1 R * 
10 C I R * 
• 0 C IR * 
7© CIR ‘ 
10 CIR * 
10 CIR * 
1 OO CIR * 


ITERATION 110 CIR 
ITERATION 120 CIR 
ITERATION 130 CIR 
ITERATION 140 CIR 
ITERATION 110 CIR 
ITERATION 110 CIR 
ITERATION 170 C1M 
ITERATION 140 CIR 
ITERATION 1*0 CIR 
ITERATION 200 CIR 
ITERATION 210 CIR 
ITERATION 220 CIR 
ITERATION 230 CIR 
ITERATION 240 CIR 
ITERATION 210 CIR 
ITERATION 290 C IR 
ITERATION 270 CIR 
ITERATION 210 CIR 
ITERATION 290 CIR 
ITERATION 300 CIR 
ITERATION 310 CIR 
ITERATION 320 CIR 
ITERATION 330 CIR 
ITERATION 340 CIR 
ITERATION 390 CIR 
ITERATION 340 CIR 
ITERATION 370 CIR 
ITERATION 310 CIR 
ITERATION 390 CIR 
ITERATION 400 CIR 


- 0 . 00402 

041317 
041144 
041177 
O 4 1710 
O 4111S 
O 41110 
0 41111 
O 41101 
041134 
0 411 II 
0.4 7002 
0.4 7040 
0 . 47010 
0.47122 
0 . 4 7 1 IS 
0 472 10 
O 47211 
O 4730 1 
0 4 734 7 
O 47393 
0 *7431 
047411 
O . 47131 
0 47177 
0 . 47122 
0.47117 
0.47712 
0 . 47711 
0 47100 
0 47143 
O 47111 
O 47*27 
047991 
O . 41001 
0 . 44047 
0 41011 
0.41124 
0 44111 
O 4IIII 
O 41234 


LOWER SLOPE 
0 11014 

0 31134 
0 19112 

0.11111 
0 . 13*71 
0 OHIO 
0 01321 
O 04412 
0 02*71 
0 O 10 1 1 
0 0014 2 
001704 
002111 
003103 
0.04731 
0 . 01494 
0.01249 
09174 
07411 
0*021 
01111 
09020 
0*413 

09 49 1 
10300 

ions 
110 12 
1 1333 
1 1114 
1 20* 1 
12449 
1 2130 
13412 


opm 

DPM 

OPM 

DPM 

DPM 

DPM 

OPM 

DPM 

DPM 

DPM 

DPM 

DPM 

DPM 

DPM 

DPM 

DPM 

DPM 

DPM 

DPM 

DPM 

DPM 

DPM 

DPM 

DPM 

DPM 

DPM 

DPM 

DPM 

DPM 

OPM 

DPM 

DPM 

DPM 

DPM 

DPM 

DPM 

OPM 

DPM 

DPM 

DPM 


AT 

AT 


00011122 AT 
0001IC13 AT 
00012100 AT 
00011414 AT 
00014043 AT 
00001411 AT 
OOO 11317 
OOOOIIII 
OOOOS734 
OOO 19192 
000167*2 AT 
00001711 AT 
00001774 AT 
00009107 AT 
.00007431 AT 
00001515 AT 
.00001213 AT 
.00001112 AT 
00021211 AT 
.00001394 AT 
00011921 AT 
00011193 AT 
00007121 AT 
OOOOIIII AT 
00009793 AT 
OOO 10111 AT 
0001027* AT 
00007192 AT 
00001314 AT 
00601717 AT 
00012022 AT 
OO00S227 AT 
. OOO 10141 AT 
.00001214 AT 
OOO 12121 AT 
00001131 AT 
00009711 AT 
00007113 AT 
00012*47 AT 
00007271 AT 


31 2 
47 M 
37 2 
4 7 10 
47 12 
31 11 
41 11 
2 1 4 
41 12 
41 1 
41 1 3 
47 13 
2 1 0 
41 12 
4 1 7 
#7 1 1 
41 2 
41 12 
4 I 1 3 
2 1 1 
41 1 2 


41 13 

41 13 

4 1 13 

47 12 

47 1 

4 1 11 

47 1 

44 12 

3 1 1 

47 11 

4 1 4 

2 1 3 

2 1 1 
41 12 

41 M 
2 10 


NS S P 
NSSP 
NS 1 P 
NSSP 
NS S P 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 
NSSP 



01 LTA Y 

OR 

YUI TE 


0 . 0391 

SEP 

A T 

X « 

0 3 7 1 


0 E L T A Y 

OR 

YU I TE 


0 . 0391 

SEP 

AT 

X < 

0 3 7 1 


OELTAY 

OR 

YUI TE 


0 . 0391 

SEP 

AT 

X * 

0.371 


DEL TAY 

OR 

YUI TE 


0 0391 

SEP 

AT 

X * 

0 37 1 


OELTAY 

OR 

YUI TE 


0.0319 

SEP 

AT 

X * 

0 3 7 1 


DEL TAY 

OR 

YUT TE 


0.0412 

SEP 

AT 

X « 

0.371 


DEI TAY 

OR 

VU I TE 


0.045 1 

SEP 

AT 

X • 

0 3 7 1 


OELTAY 

OR 

vu r Te 


0.0441 

SEP 

AT 

X * 

0 3 7 1 


OELTAY 

OR 

YUI TE 


0.04 41 

SEP 

AT 

X « 

0 3 7 1 


OELTAY 

OR 

YU I TE 


0 0441 

SEP 

AT 

X * 

0 3 7 1 


DELTAY 

OR 

YU I TE 


0 0444 

SEP 

AT 

1 * 

0 37 1 


OELTAY 

OR 

YU I TE 


0 0442 

SEP 

AT 

X « 

0.371 


DELTAY 

OR 

YUI TE 


0.044 1 

SEP 

AT 

X ■ 

0 3 7 1 


OELTAY 

DR 

YUITE 


0.0440 

SEP 

AT 

X » 

0 3 7 1 


OELTAY 

OR 

YUI TE 


0 043* 

SEP 

AT 

X • 

0 3 7 1 


OELTAY 

OR 

YUITE 


0 0431 

SEP 

AT 

X ■ 

0 3 7 1 


OELTAY 

OR 

YUI Tl 


0.0437 

SEP 

AT 

X ■ 

0 3 7 1 


DELTAY 

OR 

YUITE 


0 . 04 31 

SEP 

AT 

X • 

0 3 7 1 


OELTAY 

OR 

YUITE 


0 0435 

SEP 

AT 

X • 

0 37 1 


DELTAY 

OR 

YUITE 


0 04 3 4 

SEP 

AT 

X * 

0 3 7 1 


OELTAY 

OR 

YUITE 


0 0434 

SEP 

AT 

X * 

0 3 7 1 


DELTAY 

DR 

YUITE 


0.0433 

SEP 

AT 

X 1 

0 3 7 1 


DELTAY 

OR 

YUITE 


0.0432 

SEP 

AT 

X • 

0 3 7 1 


DELTAY 

OR 

YUITE 


0 0432 

SEP 

AT 

X « 

0 37 1 


OELTAY 

OR 

YUITE 


0 043 1 

SEP 

AT 

X * 

0 3 7 1 


DELTAY 

OR 

YUITE 


0.0431 

SEP 

AT 

X > 

0.371 


DELTAY 

DR 

YUITE 


0 0430 

SEP 

AT 

X * 

0.371 


DELTAY 

OR 

YUITE 


O 0430 

SEP 

AT 

X * 

0 3 7 1 


DELt AY 

OR 

YUITE 


0 0429 

SEP 

AT 

X * 

0 3 7 1 


OELTAY 

OR 

YUITE 


0 0429 

SEP 

AT 

X * 

0 37 1 


OELTAY 

OR 

YUITE 


0 0 4 24 

SEP 

AT 

X ‘ 

0 3 7 1 


OELTAY 

OR 

YUITE 


0 04 21 

SEP 

AT 

X • 

0.371 


DELTAY 

OR 

YUITE 


0 0427 

SEP 

AT 

X • 

0 37 1 


DELTAY 

OR 

YUITE 


O . 0427 

SEP 

AT 

X • 

0 3 7 1 


DELTAY 

OR 

YUITE 


0 . 0427 

SEP 

AT 

X » 

0 3 7 1 


DEL t AY 

OR 

YUI TE 


0 0421 

SEP 

AT 

X * 

0.371 


DELTAY 

OR 

YU i te 


0 0421 

SEP 

AT 

X « 

0.371 


DELTAY 

OR 

YUITE 


0 0421 

SEP 

AT 

X * 

0.371 


DELTAY 

OR 

YUITE 


0 0421 

seP 

AT 

X * 

0.371 


DELTAY 

OR 

VU1 TE 


0 0425 

SEP 

AT 

X > 

0.371 


BOUNDARY LAYER ANALYSIS POP REYNOLDS NUMBER OF 0 . BOOE+07 


x 

YUOR 1 C 

DU 

11 U 

YLORIC 
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